
#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkPolyDataMapper.h>
#include <vtkProperty.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkMath.h>
#include <vtkPolyDataReader.h>
#include <vtkSmoothPolyDataFilter.h>
#include "vtkScalarBarActor.h"
#include <vtkLookupTable.h>
#include <vtkFloatArray.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkWindowToImageFilter.h>
#include <vtkJPEGWriter.h>
#include <sstream>
#include <vtkCenterOfMass.h>
#include <vtkAssembly.h>
#include <vtkSphereSource.h>
#include "vtkCellArray.h"
#include "vtkTransform.h"
#include "vtkTransformFilter.h"
#include "vtkTransformPolyDataFilter.h"
#include "vtkRotationalExtrusionFilter.h"
#include "vtkLODActor.h"
#include "vtkPropAssembly.h"
#include <vtkCamera.h>
#include <vtkSphereSource.h>
#include <vtkGraphicsFactory.h> // Offscreen Rendering
#include <vtkImagingFactory.h> // Offscreen Rendering
#include "vtkAnnotatedCubeActor.h"
#include "vtkAxesActor.h"
#include "vtkTextProperty.h"
#include "vtkCaptionActor2D.h"


#define VTK_CREATE(type, var) \
  vtkSmartPointer<type> var = vtkSmartPointer<type>::New()

#define USAGE "Usage: CurvatureMapVisualization 1 2 3 4\n"\
			  "where\n"\
			  "\t1. CurvatureMap1 - Generated CurvatureMap description of the object1 (Diastole) \n"\
			  "\t2. CurvatureMap2 - Generated CurvatureMap description of the object2 (Systole)\n"\
			  "\t3. width - specify the window width in pixels \n"\
			  "\t4. height - specify the window height in pixels'\n"\
              "\t5. inSliceSpacing - Optional: inSliceSpacing = X-spacing = Y-spacing'\n"\
              "\t6. x1 position - Optional: specify the x pixel-position of reference point on object1'\n"\
              "\t7. y1 position - Optional: specify the y pixel-position of reference point on object1'\n"\
              "\t8. x2 position - Optional: specify the x pixel-position of reference point on object2'\n"\
              "\t9. y2 position - Optional: specify the y pixel-position of reference point on object2'\n"
// mkdir JPG
// ./CurvatureMapVisualization SPHARM1.vtk  SPHARM10.vtk 1024 512
void AttributeAdjustment(vtkPolyData* inputPolyData)
{

	int nTuples = vtkDoubleArray::SafeDownCast(	inputPolyData->GetPointData()->GetScalars() )->GetNumberOfTuples();
	std::cout << "nTuples: " << nTuples << std::endl;
	double* scalars = new double[nTuples];
	double mScalars = 0;
	for (int i = 0; i < nTuples; i++) // Calculate:  mean(x(:))
	{
		scalars[i] = vtkDoubleArray::SafeDownCast( inputPolyData->GetPointData()->GetScalars() )->GetValue(i);
		mScalars += scalars[i];
//std::cout << "scalars[i]: " << scalars[i] << std::endl;
	}
	mScalars /= nTuples;
	std::cout << "1--mScalars: " << mScalars << std::endl;

	for (int i = 0; i < nTuples; i++) // Update:  x = x-mean(x(:));
	{
		scalars[i] = scalars[i] - mScalars;
	}

	mScalars = 0;
	for (int i = 0; i < nTuples; i++) //ReCalculate:  mean(x(:))
	{
		mScalars += scalars[i];
	}
	mScalars /= nTuples;
	std::cout << "2--mScalars: " << mScalars << std::endl;

	double madScalars = 0;
	for (int i = 0; i < nTuples; i++) // Calculate:  mad(X(:)) = MEAN(ABS(X-MEAN(X))
	{
		madScalars += fabs(scalars[i] - mScalars);
	}
	madScalars /= nTuples;
	std::cout << "madScalars: " << madScalars << std::endl;

	for (int i = 0; i < nTuples; i++) // Update:  x = x/(2*mad(x(:)));
	{
		scalars[i] = scalars[i] / (2 * madScalars);
		if (scalars[i] > 1)
			scalars[i] = 1;
		if (scalars[i] < -1)
			scalars[i] = -1;
	}

	for (int i = 0; i < nTuples; i++) // Calculate:  mean(x(:))
	{
		vtkDoubleArray::SafeDownCast( inputPolyData->GetPointData()->GetScalars() )->SetValue(i,scalars[i]);
	}

}

int main(int argc, char *argv[])
{

	float inSliceSpacing = 0;
	int x1Position = 0, y1Position = 0, x2Position = 0, y2Position = 0, allPositive = 0;

	if (argc > 10)
	{
		inSliceSpacing = atof(argv[5]);
		x1Position = atoi(argv[6]);
		y1Position = atoi(argv[7]);
		x2Position = atoi(argv[8]);
		y2Position = atoi(argv[9]);
		allPositive = atoi(argv[10]);
	}


	  // Setup offscreen rendering
	vtkSmartPointer < vtkGraphicsFactory > graphics_factory = vtkSmartPointer< vtkGraphicsFactory > ::New();
	graphics_factory->SetOffScreenOnlyMode(1);
	graphics_factory->SetUseMesaClasses(1);

	vtkSmartPointer < vtkImagingFactory > imaging_factory = vtkSmartPointer< vtkImagingFactory > ::New();
	imaging_factory->SetUseMesaClasses(1);

	vtkPolyDataReader *reader1 = vtkPolyDataReader::New();
	reader1->SetFileName(argv[1]);
	reader1->Update();

	vtkPolyDataReader *reader2 = vtkPolyDataReader::New();
	reader2->SetFileName(argv[2]);
	reader2->Update();

	int width, height;
	width = atoi(argv[3]);
	height = atoi(argv[4]);

	vtkSmartPointer<vtkLookupTable> lookupTable = vtkSmartPointer<vtkLookupTable>::New();
	lookupTable->SetTableRange(-1.0, 1.0);
//	lookupTable->SetTableRange(0, 0.15);
	lookupTable->SetNumberOfColors(512);
	lookupTable->SetNumberOfTableValues(512);
	lookupTable->Build();
    for (int i = 0; i < 512; i++)
    {
    	double 	rgba[4];
    	lookupTable->GetTableValue(i, rgba);
        lookupTable->SetTableValue(i, rgba[2], rgba[1], rgba[0], rgba[3]);
    }


    //////////////////
	 // Object 1
	vtkSmartPointer < vtkPolyData > inputPolyData1 = reader1->GetOutput();
	if (allPositive == 0)
	AttributeAdjustment(inputPolyData1);


	// Top slice location
	double bounds[6];
	inputPolyData1->GetBounds(bounds);
	double polyData1Zmax = bounds[5];


	// Compute the center of mass
	vtkSmartPointer < vtkCenterOfMass > centerOfMassFilter = vtkSmartPointer< vtkCenterOfMass > ::New();

#if VTK_MAJOR_VERSION <= 5
	centerOfMassFilter->SetInput(inputPolyData1);
#else
	centerOfMassFilter->SetInputData(inputPolyData1);
#endif
	centerOfMassFilter->SetUseScalarsAsWeights(false);
	centerOfMassFilter->Update();

	double center1[3];
	centerOfMassFilter->GetCenter(center1);
	std::cout << "Mass center of Object1 is: " << center1[0] << " " << center1[1]<< " " << center1[2] << std::endl;


	// Create input mapper
	vtkSmartPointer < vtkPolyDataMapper > inputMapper1 = vtkSmartPointer< vtkPolyDataMapper > ::New();
#if VTK_MAJOR_VERSION <= 5
 inputMapper1->SetInput(inputPolyData1);
#else
 inputMapper1->SetInputData(inputPolyData1);
#endif

	if (allPositive == 1)
	{
	inputMapper1->SetScalarRange(0, 0.03);
	}
	else
	{
	inputMapper1->SetScalarRange(-1, 1);
	}
   inputMapper1->SetLookupTable(lookupTable);

	// Create input actor
	vtkSmartPointer < vtkActor > inputActor1 = vtkSmartPointer < vtkActor > ::New();
	inputActor1->SetMapper(inputMapper1);
	inputActor1->SetOrigin(center1[0],center1[1],center1[2]);
	inputActor1->SetPosition(-center1[0],-center1[1],-center1[2]);
//	inputActor1->RotateX(-60);


	// Object 2
	vtkSmartPointer < vtkPolyData > inputPolyData2 = reader2->GetOutput();
	if (allPositive == 0)
	AttributeAdjustment(inputPolyData2);

	// Top slice location
	inputPolyData2->GetBounds(bounds);
	double polyData2Zmax = bounds[5];

	// Compute the center of mass
#if VTK_MAJOR_VERSION <= 5
	centerOfMassFilter->SetInput(inputPolyData2);
#else
	centerOfMassFilter->SetInputData(inputPolyData2);
#endif
	centerOfMassFilter->SetUseScalarsAsWeights(false);
	centerOfMassFilter->Update();

	double center2[3];
	centerOfMassFilter->GetCenter(center2);
	std::cout << "Mass center of Object2 is: " << center2[0] << " "	<< center2[1] << " " << center2[2] << std::endl;


	// Create input mapper
	vtkSmartPointer < vtkPolyDataMapper > inputMapper2 = vtkSmartPointer< vtkPolyDataMapper > ::New();
	#if VTK_MAJOR_VERSION <= 5
	  inputMapper2->SetInput(inputPolyData2);
	#else
	  inputMapper2->SetInputData(inputPolyData2);
	#endif


	if (allPositive == 1)
	{
	inputMapper2->SetScalarRange(0, 0.03);
	}
	else
	{
	inputMapper2->SetScalarRange(-1, 1);
	}
	inputMapper2->SetLookupTable(lookupTable);

	// Create input actor
	vtkSmartPointer < vtkActor > inputActor2 = vtkSmartPointer < vtkActor > ::New();
	inputActor2->SetMapper(inputMapper2);
	inputActor2->SetOrigin(center2[0], center2[1], center2[2]);
	inputActor2->SetPosition(-center2[0], -center2[1], -center2[2]);
//	inputActor2->RotateX(-60);


	// Coordinate system
		vtkPolyData *pd2 = vtkPolyData::New();
		vtkCellArray *ca2 = vtkCellArray::New();
		vtkPoints *pts2 = vtkPoints::New();
		pts2->InsertNextPoint(0, 1, 0);
		pts2->InsertNextPoint(8, 1, 0);
		pts2->InsertNextPoint(8, 2, 0);
		pts2->InsertNextPoint(10, 0.01, 0);
		ca2->InsertNextCell(4);
		ca2->InsertCellPoint(0);
		ca2->InsertCellPoint(1);
		ca2->InsertCellPoint(2);
		ca2->InsertCellPoint(3);
		pd2->SetPoints(pts2);
		pd2->SetLines(ca2);

		vtkTransform *arrowT2 = vtkTransform::New();
		arrowT2->Scale(1, 0.6, 1);
		arrowT2->RotateY(90);

		vtkTransformPolyDataFilter *arrowTF2 = vtkTransformPolyDataFilter::New();
		arrowTF2->SetInput(pd2);
		arrowTF2->SetTransform(arrowT2);

		vtkRotationalExtrusionFilter *arrowREF = vtkRotationalExtrusionFilter::New();
		arrowREF->SetInput(arrowTF2->GetOutput());
		arrowREF->CappingOff();
		arrowREF->SetResolution(30);

		vtkPolyDataMapper *spikeMapper = vtkPolyDataMapper::New();
		spikeMapper->SetInput(arrowREF->GetOutput());

		vtkLODActor *a1Actor = vtkLODActor::New();
		a1Actor->SetMapper(spikeMapper);
		a1Actor->SetScale(1, 1, 5);
		a1Actor->RotateY(180);
		a1Actor->SetPosition(0, 0, 0);
		a1Actor->GetProperty()->SetColor(1,1, 1);
		a1Actor->GetProperty()->SetSpecularColor(1,1, 1);
		a1Actor->GetProperty()->SetSpecular(0.3);
		a1Actor->GetProperty()->SetAmbient(0.2);
		a1Actor->GetProperty()->SetDiffuse(0.8);
		a1Actor->GetProperty()->SetSpecularPower(20);

		vtkLODActor *a2Actor = vtkLODActor::New();
		a2Actor->SetMapper(spikeMapper);
		a2Actor->SetScale(1, 1, 4);
		a2Actor->RotateY(-90);
		a2Actor->SetPosition(0, 0, 0);
		a2Actor->GetProperty()->SetColor(1,1, 1);
		a2Actor->GetProperty()->SetSpecularColor(1,1, 1);
		a2Actor->GetProperty()->SetSpecular(0.3);
		a2Actor->GetProperty()->SetAmbient(0.2);
		a2Actor->GetProperty()->SetDiffuse(0.8);
		a2Actor->GetProperty()->SetSpecularPower(20);

		vtkLODActor *a3Actor = vtkLODActor::New();
		a3Actor->SetMapper(spikeMapper);
		a3Actor->SetScale(1, 1, 4);
		a3Actor->RotateZ(90);
		a3Actor->RotateY(-90);
		a3Actor->SetPosition(0, 1, 0);
		a3Actor->GetProperty()->SetColor(1,1, 1);
		a3Actor->GetProperty()->SetSpecularColor(1,1, 1);
		a3Actor->GetProperty()->SetSpecular(0.3);
		a3Actor->GetProperty()->SetAmbient(0.2);
		a3Actor->GetProperty()->SetDiffuse(0.8);
		a3Actor->GetProperty()->SetSpecularPower(20);

		// reference ball
		vtkSphereSource *ss = vtkSphereSource::New();
		ss->SetRadius(1);
		vtkPolyDataMapper *rpMapper = vtkPolyDataMapper::New();
		rpMapper->SetInput(ss->GetOutput());

		vtkLODActor *rpActor1 = vtkLODActor::New();
		rpActor1->SetMapper(rpMapper);
		rpActor1->SetPosition(-x1Position*inSliceSpacing-center1[0], -y1Position*inSliceSpacing-center1[1], 5+polyData1Zmax-center1[2]);
		rpActor1->GetProperty()->SetSpecularColor(1,1,0);
		rpActor1->GetProperty()->SetSpecular(0.3);
		rpActor1->GetProperty()->SetAmbient(0.2);
		rpActor1->GetProperty()->SetDiffuse(0.8);
		rpActor1->GetProperty()->SetSpecularPower(20);

		vtkLODActor *rpActor2 = vtkLODActor::New();
		rpActor2->SetMapper(rpMapper);
		rpActor2->SetPosition(-x2Position*inSliceSpacing-center2[0], -y2Position*inSliceSpacing-center2[1], 5+polyData2Zmax-center2[2]);
		rpActor2->GetProperty()->SetSpecularColor(1,1,0);
		rpActor2->GetProperty()->SetSpecular(0.3);
		rpActor2->GetProperty()->SetAmbient(0.2);
		rpActor2->GetProperty()->SetDiffuse(0.8);
		rpActor2->GetProperty()->SetSpecularPower(20);


		// combine orientation markers into one with an assembly
		vtkSmartPointer < vtkAssembly > assembly1 = vtkSmartPointer	< vtkAssembly > ::New();
//		assembly1->AddPart(a1Actor);
//		assembly1->AddPart(a2Actor);
//		assembly1->AddPart(a3Actor);
		if (argc > 9) assembly1->AddPart(rpActor1);
		assembly1->AddPart(inputActor1);

		vtkSmartPointer < vtkAssembly > assembly2 = vtkSmartPointer	< vtkAssembly > ::New();
//		assembly2->AddPart(a1Actor);
//		assembly2->AddPart(a2Actor);
//		assembly2->AddPart(a3Actor);
		if (argc > 9) assembly2->AddPart(rpActor2);
		assembly2->AddPart(inputActor2);


	// Create ScalarBar
	VTK_CREATE(vtkScalarBarActor, scalarBar);
	scalarBar->SetLookupTable(lookupTable);
	scalarBar->SetTitle("Curvature");
	scalarBar->GetPositionCoordinate()->SetCoordinateSystemToNormalizedViewport();
	scalarBar->GetPositionCoordinate()->SetValue(0.1, 0.05);
	scalarBar->SetOrientationToVertical();
	scalarBar->SetWidth(0.8);
	scalarBar->SetHeight(0.9);
	scalarBar->SetLabelFormat("%-#6.3f");


	// Construct orientation markers
	VTK_CREATE(vtkAnnotatedCubeActor, cube);
	// anatomic labelling
	cube->SetXPlusFaceText("L");
	cube->SetXMinusFaceText("R");
	cube->SetXFaceTextRotation(-90);
	cube->SetYPlusFaceText("H");
	cube->SetYMinusFaceText("F");
	cube->SetZPlusFaceText("A");
	cube->SetZMinusFaceText("P");
	cube->SetZFaceTextRotation(90);
//	cube->SetOrigin(0,0,0);
//	cube->SetPosition(0,0,0);

	// change the vector text colors
	vtkProperty* property = cube->GetXPlusFaceProperty();
	property->SetColor(0, 0, 1);
	property->SetInterpolationToFlat();
	property = cube->GetXMinusFaceProperty();
	property->SetColor(0, 0, 1);
	property->SetInterpolationToFlat();
	property = cube->GetYPlusFaceProperty();
	property->SetColor(0, 1, 0);
	property->SetInterpolationToFlat();
	property = cube->GetYMinusFaceProperty();
	property->SetColor(0, 1, 0);
	property->SetInterpolationToFlat();
	property = cube->GetZPlusFaceProperty();
	property->SetColor(1, 0, 0);
	property->SetInterpolationToFlat();
	property = cube->GetZMinusFaceProperty();
	property->SetColor(1, 0, 0);
	property->SetInterpolationToFlat();

	VTK_CREATE(vtkAxesActor, axes2);
	// simulate a left-handed coordinate system
	VTK_CREATE(vtkTransform, transform);
	transform->Identity();
//	  transform->RotateY(180);
	axes2->SetShaftTypeToCylinder();
	axes2->SetUserTransform(transform);
	axes2->SetXAxisLabelText("x");
	axes2->SetYAxisLabelText("y");
	axes2->SetZAxisLabelText("z");

	axes2->SetTotalLength(1.5, 1.5, 1.5);
	axes2->SetCylinderRadius(0.500 * axes2->GetCylinderRadius());
	axes2->SetConeRadius(1.025 * axes2->GetConeRadius());
	axes2->SetSphereRadius(1.500 * axes2->GetSphereRadius());

	VTK_CREATE(vtkTextProperty, tprop);
	tprop = axes2->GetXAxisCaptionActor2D()->GetCaptionTextProperty();
	tprop->ItalicOn();
	tprop->ShadowOn();
	tprop->SetFontFamilyToTimes();

	axes2->GetYAxisCaptionActor2D()->GetCaptionTextProperty()->ShallowCopy(tprop);
	axes2->GetZAxisCaptionActor2D()->GetCaptionTextProperty()->ShallowCopy(tprop);

	// combine orientation markers into one with an assembly
	VTK_CREATE(vtkAssembly, marker);
	marker->AddPart(cube);
	marker->AddPart(axes2);


	// Create renderers
	VTK_CREATE(vtkRenderer, renderer1);
	renderer1->AddActor(assembly1);
	renderer1->SetViewport(0.45, 0, 0.9, 1.0);
	renderer1->SetBackground(.1, .2, .3);

	VTK_CREATE(vtkRenderer, renderer2);
	renderer2->AddActor(assembly2);
	renderer2->SetViewport(0, 0, 0.45, 1.0);
	renderer2->SetBackground(.1, .2, .3);

	VTK_CREATE(vtkRenderer, renderer3);
	renderer3->AddActor(marker);
	renderer3->SetViewport(0.35, 0, 0.55, 0.4);
	renderer3->SetBackground(.1, .2, .3);

	VTK_CREATE(vtkRenderer, renderer4);
	renderer4->AddActor2D(scalarBar);
	renderer4->SetViewport(0.9, 0, 1.0, 1.0);
	renderer4->SetBackground(0.3, 0.3, 0.3);


    // Reposition Camera
	vtkCamera* cam1 = (renderer1->GetActiveCamera());
//	cam1->Azimuth(50);
	cam1->Elevation(-60);
	renderer1->ResetCamera();

	vtkCamera* cam2 = (renderer2->GetActiveCamera());
//	cam1->Azimuth(50);
	cam2->Elevation(-60);
	renderer2->ResetCamera();

	vtkCamera* cam3 = (renderer3->GetActiveCamera());
//	cam1->Azimuth(50);
	cam3->Elevation(30);
//	cam3->Zoom(1.5);
	renderer3->ResetCamera();


 for (unsigned int num = 0; num <= 72; num++)
 {

//	 marker->RotateZ(5);
	 assembly1->RotateZ(5);
	 assembly2->RotateZ(5);

	 cam3->Azimuth(-5);
//	 cam1->SetViewUp(0,1,0);
	 renderer3->ResetCamera();



	//Add renderer to renderwindow and render
	vtkSmartPointer < vtkRenderWindow > renderWindow = vtkSmartPointer< vtkRenderWindow > ::New();
	renderWindow->SetOffScreenRendering( 1 );
	renderWindow->AddRenderer(renderer1);
	renderWindow->AddRenderer(renderer2);
	renderWindow->AddRenderer(renderer3);
	renderWindow->AddRenderer(renderer4);
	renderWindow->SetSize(width, height);
	renderWindow->Render();


	vtkSmartPointer<vtkWindowToImageFilter> windowToImageFilter = vtkSmartPointer<vtkWindowToImageFilter>::New();
	windowToImageFilter->SetInput(renderWindow);
	windowToImageFilter->Update();

	vtkSmartPointer<vtkJPEGWriter> imageWriter = vtkSmartPointer<vtkJPEGWriter>::New();
		std::string BaseName = "./JPG/Object.jpg"; // output contours sliced from original object
		std::stringstream out;
		out << setfill('0') << setw(2) << num;
		BaseName.insert(6, out.str()); //13
		imageWriter->SetFileName(BaseName.c_str());
		imageWriter->SetInputConnection(windowToImageFilter->GetOutputPort());

		if (num)
			imageWriter->Write();

 }
	return EXIT_SUCCESS;
}
